{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1: MAX-CUT with the CIM\n",
    "\n",
    "This notebook serves as an introduction to utilizing the coherent Ising machine for MAX-CUT problems. The key goals are as follows:\n",
    "\n",
    "* Brute force solving a simple MAX-CUT problem\n",
    "* Solving a MAX-CUT problem with our CIM simulator\n",
    "* Showcasing spin and energy evolutions of the CIM simulator.\n",
    "\n",
    "The physical realization of the coherent Ising machine is a computing system that consists of numerous coupled optical parametric oscillators with increasing gain that aims to find the minimum energy of some Ising problem. \n",
    "\n",
    "Due to physical constraints, we simulate our CIM in Python to model its performance. This means that we attempt to compute the minimum of the energy function \n",
    "\n",
    "$$ H = -\\sum_{1\\leq i \\leq j \\leq N} J_{ij}\\sigma_i \\sigma_j - \\sum_{1 \\leq i \\leq N} h_i \\sigma_i$$\n",
    "\n",
    "where we have an $N \\times N$ matrix for $J$ (which is our spin coupling matrix) and an $N$-dimensional vector $h$ (which serves as our external field terms), along with each spin $\\sigma_i \\in \\{ -1, 1\\}$.  This means to leverage the CIM on some optimization problem requires converting it into an Ising problem where the solution to the problem is correlated with the ground state energy of the problem in Ising form. \n",
    "\n",
    "We provide three implementations of the CIM: chaotic amplitude control (CAC), amplitude-heterogeneity correction (AHC), and amplitude-heterogeneity correction with external field terms (AHC). The first two solvers do not have an external field, while the latter requires an external field $h$.\n",
    "\n",
    "When the MAX-CUT problem is reformatted into an Ising problem, there is no external field, so either the CAC solver or the AHC no external field solver can be used. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading J matrix\n",
    "\n",
    "The coupling matrices and bias vectors are stored as Numpy arrays, which could be loaded in from delimited files via the helper function load_matrix_from_rudy(), or from npz files.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initialization and Imports\n",
    "import sys\n",
    "from pathlib import Path\n",
    "sys.path.append(str(Path.cwd()) + \"\\\\..\\\\\") # append repo to list of search directories \n",
    "\n",
    "from cim_optimizer.solve_Ising import *\n",
    "from cim_optimizer.CIM_helper import brute_force\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "instance_path_str_MAXCUT = str(Path.cwd()) + \"\\\\..\\\\instances\\\\MC_Instances_NPZ\\\\\"\n",
    "\n",
    "# 20 spin MAXCUT problem\n",
    "N = 20\n",
    "mc_id = 1 # select first example of 20 spin MAXCUT problem\n",
    "J = - np.load(instance_path_str_MAXCUT + f\"MC50_N={N}_{mc_id}.npz\") # load J matrix for 50% density MAX-CUT problem\n",
    "gamma = 0.010 # set gamma hyperparameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we only modify the gamma hyperparameter since the default value for gamma (which is meant for MIRP problems) is too low for convergence towards the ground state energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brute force solving"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We test every single possible spin configuration $\\big($ i.e. $2^{20}$ configurations since each spin $\\sigma_i \\in \\{ -1, 1\\} \\big)$, and return a tuple containing an array of the minimum spin configuration and ground state energy (lowest energy achieved)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The spin configuration in the ground state is [-1 ... -1]\n",
      "The ground energy is -29.0\n"
     ]
    }
   ],
   "source": [
    "spins_ground, E_ground = brute_force(J)\n",
    "\n",
    "print(\"The spin configuration in the ground state is {}\".format(spins_ground))\n",
    "print(\"The ground energy is {}\".format(E_ground))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving with the coherent Ising machine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now apply the CAC CIM model to our J matrix for a single trial with default hyperparameters (with the exception of gamma). Our code returns certain pertinent metadata, including the target Ising energy, best Ising energy found and the corresponding spin configuration, time elapsed, number of trials, and success probability across the trials.\n",
    "\n",
    "Note that the default behaviour of `Ising().solve` is to attempt to optimize the hyperparameters initially, but we showcase the highly efficient single-run performance of the simulated coherent Ising machine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No External Field Detected\n",
      "Target Ising Energy: -inf.\n",
      "Best Ising Energy Found: -29.0.\n",
      "Corresponding Spin Configuration: [1. ... 1.].\n",
      "Time Elapsed: 0.32407522201538086.\n",
      "Number of Runs Completed: 1.\n"
     ]
    }
   ],
   "source": [
    "test = Ising(J).solve(cac_gamma=gamma, hyperparameters_randomtune = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metadata and internal dynamics of the CIM\n",
    "\n",
    "Various properties of the solver runs can be read from the attributes of the run result object, such as the time elapsed, best Ising energy reached by the CIM, as  well as trajectories of the spin amplitudes over time. Similarly, the evolution of the Ising energy corresponding to the spin amplitudes over time could be plotted by specifying the argument plot_type=\"energy\". \n",
    "\n",
    "If the solve() function was called with multiple runs, and specific runs among them need to be plotted, the argument trajectories_to_plot=[] takes in a list of the desired run indices.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time Elapsed: 0.32407522201538086 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Time Elapsed: {} seconds\".format(test.result.time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum Energy Achieved: -29.0\n"
     ]
    }
   ],
   "source": [
    "print(\"Minimum Energy Achieved: {}\".format(test.result.lowest_energy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy Evolution: [[-17. ... -29.]]\n",
      "Spin Evolution: [[[2.1365203e-04 ... 9.9761528e-01]\n",
      "  ...\n",
      "  [9.9497673e-05 ... 9.9525070e-01]]]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(threshold=10, edgeitems=1)\n",
    "print(\"Energy Evolution: {}\".format(test.result.energy_evolution))\n",
    "print(\"Spin Evolution: {}\".format(test.result.spin_trajectories))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The internal dynamics of both the spin amplitudes and energy evolution are included to showcase how our CIM evolves over discrete time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAACZCAYAAABNGxy5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABHcklEQVR4nO3dd3wc1bnw8d+Zme276l2yJDfcwQ3bdLCpCTUhlIQQSCHJS3LTc28u6fWm90YSCAkkJKEFCL0YMMXgggvuli0Xyep968yc949d2XKXsVUsP98Py+7MnJl5dsYjPToz5xyltUYIIYQQQowMxlAHIIQQQgghjh1J7oQQQgghRhBJ7oQQQgghRhBJ7oQQQgghRhBJ7oQQQgghRhBJ7oQQQgghRhBrqAMYLgoKCnR1dfVQhyGEEEIIcVhLly5t1loXHmiZJHcZ1dXVLFmyZKjDEEIIIYQ4LKVU7cGWyW1ZIYQQQogRRJI7IYQQQogRRJI7IYR4G3bVdfDA35ayctmOoQ5FCCH2IsmdEEK8Dcvf2EZXZ5zlr9fS1NA11OEIIcRuktwJIcQRslMOjbu6mHxyGT6/R2rvhBDDirSWFUKII9TeHgOtKSqJYFkmK5fvoKszTiTLP9ShCSGE1NwJIcSRam+JUloaoCQ/xcSJYZSC9W/tGuqwhBACGMHJnVJqlFLqeaXUGqXUW0qpTw11TEKIkaGrM87UabkoBSYJTplezMZ1DdgpZ6hDE0KIkZvcATbwOa31ZGAecKtSavIQxySEGAHisQR5eT5WPvU4//niBXQ+9WmS0S5qNjYNdWhCCDFyn7nTWtcD9ZnPXUqptUA5sGZIAxNCHP+0Q6KrjeV//gRol+btGynMeYE1vrsZO+FSTHMk/90shBjuToifQEqpamAGsHiIQxFCjAAeCzYtvBe0S0W+STDgp7k9gXrhRta8sXqowxNCnOBGfHKnlAoD9wOf1lp37rPsFqXUEqXUkqYmuZ0ihOgfy4Kdy57BVJCbbeGf+t/4I7k0tcWovf0i2pvbhzpEIcQJbEQnd0opD+nE7h6t9QP7Ltda3661nq21nl1YWDj4AQohjju27eD1Ktp2bMJrKXyBAKGxCzjnS/djeQM0NLby6tfOwHHcoQ5VCHGCGrHJnVJKAX8C1mqtfzLU8QghRoZ4LIXPckhGu/D6PFh5Yzn/utksfTXJRd94EGVY1G/bxMKvXoLWeqjDFUKcgEZscgecAbwfmK+UejPzesdQByWEOL7FYylUohUAE4eFC5t55allTJwznmiqmNM/9RdcbbB95UJe/tktQxztHq7rsvj+3/Knj8/n3i9dw1vP3S/JpxAj1IhN7rTWi7TWSmt9stZ6eub12FDHJYQ4vsVjKVKd2wEI+DSN7X5+8oU7GH9KKdu2xKmcPpcZN/4UWys2PvcXFv34Q0MccdrCP32LJ37xBVw7RcPm1dz3jQ/w509eSOOWtUMdmhDiGDtscqfSblBKfTUzXamUmjPwoQkhxPATiyWJNWwAwGspzrn2Mpq372DLv77M3Px/s/qpF5l8/mVMu+5HOFqxYeHdPPO/5w/pM3iNW9ay6J4fM/2S9/Ph21/kk/eu5rIv/Irm2g384SNn8vLffoLrSAfMQowU/am5+w1wGnB9ZroL+PWARSSEEMOYnbTpqq8BwO9VzL5kAZ//iJeC5DKclo1M1HfT9PLfmXrJlUy94XbAZOvKl/jPR6rpaqofkpiXPHQ7psfLBR//NkopDMNg5qU38fG7XmfcvAt55vdf5c5PnE9T7fohiU8IcWz1J7mbq7W+FYgDaK3bAO+ARiWEEMOUY9t0N23HUGB5TJI6zLmz4jz5ejbfuu9clm2tJLvjGbpf/g1TzzuLce/7J75AmIaGBv7z8bGseuCXgxpvMtbDyqf+weRzriSYnb/XsnBeMdd862+86yt30LJ9M7d/6HReuffnUosnxHGuP8ldSillAhpAKVUISBt/IcQJSbsO0bZd6eTOVLTu2IyBi6p+B9FoksT4W3hh8wz8HcuIv/wDZp87jpIr76Nw7Ey6Yi5L7vw8//l/E2nasm5Q4t346hMkop3MeOcH0FoT7UnS0REnFkuhtUYpxbTzr+Hjd73O2Dnn8/Rvb+OO/zefbateHZT4hBDHXn+GH/sF8CBQpJT6DnA18OUBjUoIIYYprV3inW0YhkIF8skxW9CYjD7rHJ5//m/k5JtM+q+vsOiuO5jjPk7shW9xxtmfYV3xr2Hxo7S+/gPqarfwzOemUzz5TGZ/8i6yissHLN6apQvxBbPIrprB6tW7SCb21Mp5PCY5OX7y8oOE84q59jv3surpf/DM777MnZ+4gPFzL2Li2ZdTWD0BbzCCQuHYSVLxKKlELP3e+zkzbSfiePwBwnklhPNLKKyeRKSghHTvVEKIwaD60xReKTURWAAo4Fmt9YhrXjV79my9ZMmSoQ5DCDHMrXx9PUu/OxuP4TJpximMP/0CWmo2sLXq6/gCFnd+6y5+/sCX8Af9LH7gUSo7/ojfZ+Cd/XFivnE8/tc38NX/lq5tr2EoTSSgKBg3l6k3/YKiCacc83h/cd00skpPYtr7fk68J8bz9y1k44oaSqpLmD1/BlUTqzAtA5/PQ15+kLy8AIZO8to/f8WSf/+RrpbDPCfYj95UvOF8ssqnEi6bBVmn0N4GF757JpNmjD42X1KIE5BSaqnWevYBlx0suVNK5R1qo1rr1mMQ27AhyZ0Qoj9WL17Jq9+cQ3bIw9h5FzJ+TA5vbYDE1M8SLojQXN/Ki/c/zzf+9Ek8HouaZSvhje+SF2rHKVtA4OSrqVnXyXP3PENOx+9JdGzDVJqwH3KKRlF89i1MuvJTeP2+o461YetGfnfjDKZcdRutqSn88Rt3MW5qFdPmnESsJ87aZZvZUdPAjHNO4bx3n0315CoMQxEKecnLD5GT7aOzYSttO2tIJWPYKRvHNXHwkEiZxBMK0xvAsPwYXj9dbTZdbSkSdW+iWpZiRDdiR+uIRdvp6IoRi9soQ1E5fjJVF32Z9TUBZp05hqkzqwgE049yR9saWPPo76hf9SJ2PEqosIKiiXOpnPMO8qomH/UxEWKkeLvJ3RbSf5MpoBJoy3zOAbZprUfUn1yS3AkhDsd1NWtefI5Xf/hOinIsqs+9iergFl5YdxLTb7qNgpIsVq+qo7WxnRfvX8gXfnwzuQXZdDa3sf4f32N81kpSniL8Mz6Ap2giG1Y18eyfHyTc+XeIbcNQmqAXIgGTUMXJlJx5M+MuvBGPL3BEcWqtaW+L8dI/72D53/6HBt9N1G5J8OnvfZDuNsX2jU2kkg65hWH8IZPGul2seH01tRt3cOr5Mznz0tMoqU7fSjUUoMB192wbDS272thV00LT9i52bWnB07GSU6rWUJ27GS9tKMAMZKPDlUTdbKJJL9FolIYtK6jb1YrPA+PLAgTCudiuhfKGUR4fzdvXYydj5FdPJVhQQVfDVtq3p1vx5lZOYsxZVzPuvOuIFFcf03MrxPHmbSV3fVb+A/BgbwfASqlLgCu11h895pEOIUnuhBCHE40m2fjM3bz++//HqAKT0Zd8ilJnOfcvm8uN3/8qpmUQjSZZtWInqYTNY3c9wZgJxVz1wQvJzc/itX/eR3nXveQEo6jiU7EmXY0RLqS1oYfn7n6Yljf/TtB9C9D4rXSiF/Qb+PJHkzNpPlULbiZvzAyUceC2cFpruruTbN/WSrTH5rU7/4ddKx6Bsd8jGChi1Wu1aFdTPCoHj9eirbGbnq747vVDWX5CWRZdXe3Ek1EmnlpNdkEEj9dLtDNOT3uC9l1xWnZGScZSlJS4zJ+zjVHeV6BnG6Y3TGTCuXiqzqU+OZ7NOzwkUw7ZOUEq8jtJvPZ/dKx/ge6URU2TSzLlML5YUVw+ih4nQuO2dSjtkBvSeE0DZfnIGj2bnLGziKVcdq1bSsPa1wAomXom4+e/l9FnXIU3mDUYp1+IYeVok7tVWutph5t3vJPkTghxOC1N3Wx65LusvP8njC31MPayz5PT9Sr3rb2KD31nz9+7yaTNmlV1pGxoqmtm0cOv0NnSTm5+mPzCHKq9qzlj3Fo8lkt3aDrW+HeSUz4G0NSsWMcb//oT8dqnMHQ3hmEQsFyCXo3HBJ/fh5lTia9sFoGquZh544mrQuKOhenz4fX56Grv5um/LURv/hloi60dN5NbFObcK6dy7lXTKB6VszvW7vYY2zc1U7u+iW0bm6hd18SOTc0k4inSvx56f0coPF6TiTNLmTu9gzLjBZJbn0W7NuHRp5E35wZ6cs5i/YYO6ne0YxgGVWPyGDcui87Fv2XLYz/H9AYYe9lnKTvt3cTjMR789sfYvmYZY0stAqZLe4/mnM/8Ee0vp3PnWlrXLya+YxmqcyMKjeHxE6qcQcoM0LRtE92N2zG9fqrmXcr4895H+Yz5GGZ/2gkKcfw72uTuSeAl4O7MrPcBZ2utLzqmUQ4xSe6EEIezfWsLm//xSTa++ACTq7yMu/zz+Jtf4eGmT/P+z1y8V1mtNR3tcbZuaca209O2bdO4vYl4NEHQjDPaWky5fh0Dm2bGkyg4g+xxs8gvCKKdFDvefJENLzzArhXP4TopvD4fAdPGZySxDDBNhddr4vOZeML5pIw8uqMmne0Qi1tsqV9BYeEo5p17DoWjCjB9fsxADla4ECtcgCecjxnKxwrlYXhDu1u0uo5Lw/Z2tm1opqOlh3CWl+LsBnzti+hY8QCJ1m1YoXzyT72e4LRr2daSxca1DfR0JwiGfZw0qZhxE4voWvMEa+7+IvGW7ZSedjVjL/8cPq8JndvRPbuwO+t47sGHWL5kNV5TM6HMYHSJl1BuIZHiSoKl4zDDJaRUkIYdTTTV1tC1Yy26sxatNY4ZxAkU0tnajJ2IEcgpYty51zFu/nvJHz2i6h+E2M/RJnd5wNeAszOzXgS+IQ0qhBAnmk3r69nwx2upX/s608b6GXfJrXRvX8Nb2V/lkmtmHHS9ZNKhqytOPGYTiyVJJmxs2yHanaSnuZGSxEuUuq9jud3YVi7d2XNxCmbgzasiFDJRqQ62LX6cmpf/TdPGZWit8UZK8FmKoNOISQoAwxfB9IXw+nzY2uXl5TuYNbmAyiIv2k7h2nG0ax8wRmVYmIEczEA2ViAbM5iDMj3YPa3Ed63FSXSjlEl43Jnkz30/8bwz2LC+jW1bW9GuprQihwlTShhVlUescTNr/vp5mlY8RahkLCeddxVZfhu3bTMkOzM7NFHBAoxwMQ/d/SAbN20nFI5w+mkTCPRsRqeiKGXii+QSCAXwB4P4Q2FMj4eUDW2dNq1N7XTsqiPR3UrchljKIJHSaCCUV8qoUy9m/PwbKJ4076C3soU4Xh1VcneikOROCHE461ZtZ/Uvzqe7oZap40KMOfs9bK2JoU77OjNOe3ttzDraovzrD6+wdskW5kzpYP6MBvw9a0C7qFARVvlczKJpGPkTwfLTUbeZmpfuZ8srD9OxcxOG5aFo/HSyCspINm4msWs1WmuaUvnU1rcw4do/UzVtJgVFEUxDoZNd6EQrKtmBSrRhpNox7U4Mpxs30YkTa8eJtuPEOnCdJFYgB3/xBPyjZpHKn8fOBs3Wmhai3Ql8Pg/jJhZx0qRisnICOMkYm/79fWoe/QmGoaiacBJlFXkow0CFS7GKpmIUTMTMn4iRXYkyLBJdbfztA2MpnX0pG1evoG1nDTPeeSOnnLmArk2LiG1+EdW1BdXbd77hxfRn4Q0E8fosLGXjOA7RngSxaIKuzigtrT3EkppEJo81TQN/pIBA4QR8o87EKpiKlV2OFcjGsgwMS2HHOoi11JLoaiQV78SOd2HHu7GTUcBFuzp9TgwDjy+ANxjCFwzjC4YIZmUTzi0glJNPMDuPUE4+oex8TI/n2PzDE+IAjrbm7nkO0JOR1nr+sQlveJDkTghxOGve3MKy781FJ7uZPKWasZNPYvGmCk6++VuUjso9qm2vXrKNB+9cTEtjF3NOL+Id81187cuwd60AJ93owcipTidGOaMhu4q2lnZqFj1EzaIHibU3kT9mGqdcdStGspun//Bt6ne1MmtKFdHRHyaaNQ8O05Gw5TEJBDz4/B68XhPLY4KG7u4E7a1RXNfFMA3KR+VQPbaAytH5WJaB1pr6F+9g3b1fIdbZSlFZEaOnTSVYNRerbDZW6SyMcMmBj+ljf+CV336aK3/2MpHSsbxw53dZfN9v8AbCzLriQ8x510exfEHaNr9J++YlRHdthp46jFgdRHeh7Ch7/4rSmVvgmnjCobPHpSsG0aTGdTV2Jkd0NbhakXIV8RTY+464pkApA8syQSlUegYAjuPg2AeuAd3reFomHo+FxzLTn00Dy1JYpsIyAO2gtIt2HdAuruPgahdcjUanE0oFu39N60wXFpkWzJn/oVGZBQY6czR07zKt0tO738HVGq3BcdOfew+fVnsfyt78QCmFUum9pQ+BTh8PdJ9/Uhql93zuez72/qzoswq6z3K9p8gB9JmpD7GHvsv6fN7z5Ojh9KPv38OWgIKCHG76x8COJX20yd2sPpN+4N2ArbX+4rELcehJcieEOJy3lm1i8ddPIeAzmDRzBmMqQ/xn5Uze8+1vphOho5RK2jz38CqefnAlSikuuOpkzn3nBMyuLdgNK3AaV+G2bkInuzJrKFSwAB0oZGtNPW+98hLd7S2MmXsBm7fswkkmmFSq6KlbT2TMXKqv+BrZY+fi6nSi49gu8XiKeGzPK5Z5TyUd7FQ64wmGvOQVhCgsjlBSlo3Xt6fRQtvqx1lz16dpr6slFAkz/pwryD/1OjyjTsf0hQ77nR/+4nzsWA9X/eK13c/8NWxezYt3/R9rX3wYtKZi6lxOOu1iyiefStGYKQSz89EaUimHeE83PS31dNTvINZST6y9nu6mncQ7G0h0tZDoaaWno5Xuzk56orE9iRIaU4Flgs8EnwcCHk3Ak562TDDUwfNhrXsTJXA02G46QbTdPZ9TriLl7Jnv6HSXMs4g3zDrTcr2e1d7EpXDDSDSe9x6E7DeJLLvsr7LYZ/Eap95B4uz7wfVz2X9WrfvDL0nhzzQ+5HQB4pGQ0FekI8+1HKEWzsyh0ruDtusSGu9dJ9ZLyulXj8mkQkhxHHEibXhOA4+j4U3nIvWCVK+kmOS2AF4vBYXXT2D2WeP46G7Xuexe5ex+PmNXPmBuUydfR1q6vVordHRJty2Gpy2GtzuOnRXPaNLFBUXTGLN6lrWvvYU9R0wZWIpM+adRkNdCbUrl7Lqp++gaPrFjLn8C+SOn3dUQ4K1LL2fTfd/jebaTXi8HiacdxUVl30VT24VhifQr21HW3fRuO51Zr73tr3KF4+dynu+eTetO2tY+dS9rHvpEZ79w9d3LzctL5H8Ejz+YDor0Zp4dzuxrnbsZHyvfXgDYfIrxjJ66rnkV55EfsU48ivGEMzOpW3LmzSsXUzThqW01KzATsToSWl6bEUwK49gVjZenx+P14NlmZimQrtOpibLSCfIjo1j2zi2QyKRJBmLEo9F0dFuDFfjs9idXQRyCggXlBMuHEUwvwJPVj6eUH46CTY92KkUqXgMO5XETiVxkgmcVBIMA6XM9O1tw0AZJpblwRsM4g+G8QbD+AIhPP7g7pfXH8y0HHZxEjHsRBQ7ESMV78GOR7ET6aHj7Hg3drwnPT8RxUmlMEwTwzBQpoXR+/J48fhDWP4wnkAEyx/CG4hg+YNYvmB6vj+E6QtgmPteD3vOrdYa7bqZ42igTDNzPI9+eLrebcM+7zpdE2qYFsqw0sdwhA+Hd9jkbp+RKgxgFpA9YBEJIcQwZbdtAcBnOXgDASCBJ7fqmO8nvyjCh76wgPUr63jwztf40w+eYcIp5Vx101xKKnJQoSKMUBFWxby91tN2nHk9jYT/80ce/M0PaWrqJO4ppXRsHoW5Bjs31bBj1VM0Ln+MUFEVJfOupmjWFWRVnozp9R82rljzNuqf/w07F/2dzuYGPB4Po0+7hMorvkKgeAKGJ3hEvzRrFz8KWlN92hUHXJ5XPoZzb/5fzr35f+lpb2bXxhU0bVlLd2sDXS27sJNxtOOgDAN/OAd/OJtAVi45pdXkllWTU1JFMKfgoDHllI1l9BnvBsB1bNq3r6N1yyo66jbRsWMDHTs30d7SSKyjKX3r9BAsX4BQfjmhskkUFZQTyiRxkeJqwkVVhItGYXqOftSR451SCgyD/g1tf+TbVr2J5bH5e+u41Z+ju5Q9NZY2sAX40EAGJYQQw42dckg0vgWA36vw+rzEUx5yysoHbJ8TTi7jCz+8kkVPruXxfy7nB59/iLMvmcxF75m+e7iuvpTlx8yuhPxJ6dodw+CF/zzDJd96kIBlMrZjC1X1y6lf/CC7Nq5l88M/ZtPDP0Ypg2B+KcHCarzZRXizirD8YQCceCfR+g107VhDT1sTAOHsbMafdzUV77wNX37/a+r2Vfvao2SVjiG3ctJhy4ZyChh76gLGnrrgiPfTH4ZpkVc9lbzqqfst065LoqedZHf77lonrV1Mjw9vMAtPMCKJmxhW+pPcTdJa71XPrZSSf8VCiBNKV1eCWEN6GCy/V2Hh0NwdoWTCIYfhPmqmZXDOO6cw88wxPHbvMhb+ZzUrX9/KDZ88hzETiw+4Tt26pXgDId7xzft44quXs+hXn+GC2+7F9QQw8k6icsp1VEWbSdQuomXlE3TVbyba3kZs5wo6t6RIJZO4bvppIkMpfAE/gUiYkjnnkD/zCnJmXovpC6NM79u+vZXs6aBu5UKmXHbrsL9FpgwDfyQPf2Rgz7UQx0p/krtXgJn7zHv1APOEEGLE6u6I0tO0DVOBx2Ng6R6aukqYUJEzKPuPZAe49qNnMOe88dz98xf45dce49LrZzH/imn7JUc71y6h9KTplE09i9M+8kNe/s2nWHH/T5hx3ZfQThI3FcMNFeGbdAVlk65It4qMt6GjzehUNzrRhdYOoFCeEGb2KIzsakxvAGX5j0kytn3Jk7h2iurTLjvqbQkh9nbQ5E4pVQKUAwGl1Az2PBGZBQQHITYhhBg2EvEEPc07MU0wvUF8RpyWnmwKSwd3XNPRJxXx+R9ewT9+9zKP3LOE5oYu3vOR0zGM9I9oO5WkYdMq5rzrYwBMvPhDNK5bzLK/f5fiiXMpn7EAw0rffNG7u+HQEMqH/HFk+gDZ0/WHGpiHz2tfe4RAbjFFE+Yc820LcaI7VM3dRcBNQAXwkz7zu4D/HcCYhBBi2EnGE8Q62/B7LKxIus+2uHfUMWspeyQCQS83fvpc8osjPPvQSjxek6tumotSisaat7BTCcompm+uKKU44+M/p3nTchb+9CNc9bNXCOaVZJYZKHPwR25wUgm2L32KMWe/R0aOEGIAHPSq0lrfpbU+D7hJa31en9flWusHBjFGIYQYcsmuJuxUEp/XxJeVBxp8ReOGLB7DUFz63lmce+kUXnxsDQv/k27sUbc23XtV2cQ9XZRa/iDz//uvpKKdLPzphzPdRAyduhUvkIp1Uz3v8iGNQ4iR6qDJnVLqhszHaqXUZ/d9DVJ8Qggx5FxXk6x/A4Cgx8br89IZD1BUVTmkcSmluPz9c5g2p4pH71nC9ppm6tYvI5hdQE7p3l205FZOYt5Hfkjdm8+z4r4fDVHEabWLH8ETCFN2yjlDGocQI9Wh6sN7uxYPA5EDvIQQ4oTQ3tZDd80iAMJBE48JuzpyqBiTP8SRpWvwrvvYGYQifv55+yvsXLOEsokzD/ic3IQLb2LM2Vez9J5v07Dm1SGINt2tSO3i/1Ax60LpPkSIAXLQZ+601r/PvH9j8MIRQojhp62pi5ZNy7BM8PtNgt4k9Z25TKs6uvFkj5VQxM9l75vN337xGDlNa5h6/nsOWE4pxZm3/pLmjct47kc3cdXPXx307j0aN7xBrK2BqrmXDup+hTiRHKq17C8OtaLW+r+OfThCCDH8RDva6WzcQSjgw4qUYlkWUd94/IH9OxIeKrPOGsvTf2kg2aCpPPnMg5bzBrM47wt38cgX5/PCTz7MBV/+Z2aYqsGxbfGjGKbFqNkXDdo+hTjRHOq27NLDvIQQYsSzbYfoxsdwXZdIEIIF5aQck5wxM4Y6tL0YhqK8oBGtPHQmD9y5ca/C8TM57ZYfsn3Jk7z6+8+l+7kbBNp1qVn0AKXTzsYXzhmUfQpxIjrUbdm7+k4rpbLSs3XXgEclhBDDxM7aFupe+QtKQV4Egl7N5sZiTr7mpKEObT/RXSsgOI7Xnt3MtFPHHLLspEs+QnfjNlbc9xOCeaXMuO5/Bjy+XWteoWvXVmZef9uA70uIE9lhOxhSSs1WSq0CVgKrlVIrlFKzDreeEEIc77TWNG9YTHPteiIhH4GsAsIhk52pKZRWDo/n7Xp1NO6gcesayqacyboVO4l2Jw67zuz3f4Nx513H0nu+xfJ/fH/AY9z03N/wBMJUn37FgO9LiBNZf3qPvAP4f1rraq11FXArcOfAhnVsKKUuVkqtV0ptUkoN/J+lQogRZe2bNWy89xMooDRXEy4opyMWYuz8K4c6tP1sfPVJAOZd8R4c22XV67WHXUcZBmd/6veMO+96lt79TZb89RsDdovWjkepWfQA1adficcfOvwKQoi3rT/JnaO1fql3Qmu9CLAHLqRjQyllAr8GLgEmA9crpSYPbVRCiONBKmnzxiP/ZuVPz6OruYGiXA/ZhcUU5PvZ6r2c5jabhU9tpK01OtSh7rbx1SfILa1m2lmnkVcUZmU/kjsAw7Q4+1O/Z8JFN/PmP3/Awh/djJ2MH/P4al66j1Ssi5MWvO+Yb1sIsbf+NJF6QSn1e+DvgAauBRYqpWYCaK2XDWB8R2MOsElrXQOglLoXuAJYMxTBLL7nZ2x76odHvuKBhnQ84B/WA/PX9kFHlFSaA/2Br/Zb4+Bx7V7/IDvp/2iWffeh+sxTBynTPwfdv9b7LTzk1vcvfpQB9HvPR/XPov9n8u3s9ABl9pqlDvLPXB8irgPvt/8VUQrbdXFsl1QySSrlYihNcY7B6KpsysZUss3/TuKBORTkBejuSrBqeR1nnjsGYwiG8Oor3t1BzZLnmHX5BzEMg4mnlLN0UQ2O7WJah4/NME3OvPWXRIqrWfKXr9FZX8N5X7iTrNKxxyQ+rTWr/v1L8qqnUjL1rGOyTSHEwfUnuTsl8/61febPIP3TdP4xjejYKQe295neAcztW0ApdQtwC0Bl5cD2NN/91kO8tbllQPchhDh2DKXIDZuEi0qITL+aDb75RCqncdroXLKyAzQ1drNs8TZ21XdRVpE9pLGuffFh7FSCqQvS/dudNK2MV55ez7bNTYyecOiWs72UUkx/z+fJqZjAi7/4GA9+6nRO/9hPGHfeew/YIfKR2P7GE7TVruGcz9x+1NsSQhzeYZO7zPiyI5LW+nbgdoDZs2cPaF8Aped/jrOtXx5gyQF+0GVm7fvsizrAp/1m9edbKNXP8uqgu+u17+oHrmg02L8m7SBxHXTL+yzqW7R3vQOtovaZ0Jreuh/Vv5UOvPu3/Qtq7/W01uy3KaUO+vX3m9135cMd4sOeggMd/70ONLp3/j7HVR3sdCl2r7OnVje9HUPvsyG19z4UKv31+nzHPcd+77jS1aOqz5nt/berd8/pWzS9zd6gNWgbXI3HUpikiHU00N7YQO22Zl5aUsemur9x+Q9uobiydPe5LygMEQx5qd3SOuTJ3aqn/0Fe2RjKJ58KwPippaBg/cq6fid3vapPu4yCsaew8Ccf4oWf3kLta49y2kd/TCi/7G3F5joOb/zla2SVjWXs2de8rW0IIY7MYZM7pVQOcCNQ3bf8cdCJ8U5gVJ/pisy8ITF5/mVMnn/ZUO1eCHEEbNvhrcUruDi4jKX//CkvvLqFh764gOt+9ybZOdkopVBKUTk6j3Wrd9HRHiM7JzAksTZv28iWZQs59+bbdieeoYifsqo8Nq/d9ba2GS6q5B3feYJVD/2cZX/7DnW3zmbOzd9mwoU3H/EfNuuevIO22rc474t3YVietxWPEOLI9OdBkcdIJ3arOL46MX4DGK+UGq2U8gLXAQ8PcUxCiOOAZZmcfPoMNibmMe/6z3De3DIa6hp49U9fIpnc056srCIbwzDYsa19yGJ944HfYVpeZl3+wb3mV48vZPvmFlz37d2UMEyTU979Wd71i8XkjzmZRb/6JI/ddgkddZv7vY3OXVt4489fpmz6eYw5891vKw4hxJHrT3Ln11p/Vmt9p9b6rt7XgEd2lLTWNvAJ4ElgLfBPrfVbQxuVEOJ4oZTipOkTaLfGM/3Cqxlb5mXZY3fT1tSw+5EJj8ekpCyL+p0d2LY76DF2tzXy5uN3M3XB1YTz9r79WjmukHg0SVN9x1HtI7t8HO/49mOc+Ylf0lKzkgc+OYc3//Wjw7aoTfZ08PS3r0GZFmd94lfyrJ0Qg6g/yd1flVIfUUqVKqXyel8DHtkxoLV+TGt9ktZ6rNb6O0MdjxDi+OL3e6jrKcZXNJ7pk4uwUzZrnvoLKdvZXaaiKgfHdtlVd3RJ1Nvx8t0/xk7GOfOGL+y3rGpcIQDbNjcf9X6UYTDxog/y7l8voWLmBSz5y9e47+Mz2PjcPTip5H7luxpr+c+XLqZjxwYW/PdfiRRXH3UMQoj+609r2STwQ+A29noCmUOPbSOEECNASWUZiZYCisfPIP/1ejY+fy9zrv0sXk/6x2dOboBwxMfWza2Uj8oZtBqqxpo1vPHg7Zxy8Q0UVI7fb3lxeTZev8W2jU2ceva4Y7LPUH4ZF9x2LzvffI7X7/wyL/z0Fl6/88tUzbuUgnEzMT1eGtcvYeNz96CU4oKv3Ef59OHaoYIQI1d/krvPAeO01kf/558QQhxnCosiNLUUEi4dT3WJj6UbNtHZvItgoArDSDesGDO+gJXLdtJQ30VJWdaAx+Q6Dg9//+P4wtksuOXrByxjmAbl1fns2Np6zPdfPn0+V/50ETuWP8v6J++g5qX7WffEHQCYHh+jz7yK2Td8jXDRwHYxJYQ4sP4kd5uA4dMNuxBCDCKlFC3RCPm5lZQU+GBDlLo1r1NQWo7Xm279WVKWxeYNzaxf00BBURirHx0HH41F9/yYneuW8q6v3EEot/Cg5coqc1m6qCbT5c6xrVFUhsGoWRcwatYFuI5Dd9M2tOMQLhqF6fEd030JIY5Mf34C9QBvKqV+r5T6Re9roAMTQohhw58Llp/ishJMU1G/etFez90ppZhycinxWIpN65sOuhmtNbbjkErZJJIpEokUiWSKVMrGdpx+jeu64dUnWPinbzF1wXt2d1p8MKWVucSjSdqae/r/Xd8GwzTJKhlNdvk4SeyEGAb6U3P3UObV14B2+CuEEMNJTl6YZHsET3YZhTmbqXvrZeyUs1eNWG5+kIqqXGprWigqCZOXHwLAdV2SqXRC5ziHb1FrGAaWZeCxTCzL3KvGbefaJTzwjZspGX8yl33x14etjSurSrd9q9/WRl5h+O1+fSHEcaY/I1Ts1e2JUmoU6T7jhBDihJBfGCHWmY0vu5SSXJOVWzaQjEexQ348lrm73IRJRbQ297Bi6U7mnlUNWpNMpfvFM00Dn9eDaRkYhsLYPdKGxtUa19U4jovjuKSS9u7+9DweC5/XomHTm9z9+SsJ5hRw3ff+idcfPGzcJaNyAKjb1sqUWaMOXVgIMWL0p+YOpVQh8B7geqAMeHAggxJCiOHEYxm0uiGysgvJz7Zw3TgttesIZ83ZK7mzPCanzCpn+7Y24vEkhqHweS28Xg+mebCnYFT6+RgTyAzgoLXGtl1Stk0qZbPhlcd5/Ae3EMop5P0/fZSswvJ+xR0M+cgpCFG/re1ovr4Q4jhz0GfulFIRpdQHlFJPAq8DY4HRmT7jPj9oEQohxBBTStGT9GKF88mLmKChpfYtUvbez8k5josyoGxUFi2NPdRubsPn8x4isTv4/jwek4Dfy5on7uSRb99IbvlYrv7+I5jhIlIp+/AbySgdlSvJnRAnmEPV3DWSTuq+DCzSWmul1FWDE5YQQgwvcSeI6c8mEvLi95k0b1mN62ps28GyTJJJm1g8iVKKcMhPpzdJzYYWTNNgysmlR9xaNdbZxsPf/zjrFj3KhDPeyVW3/RHlCRBPJOmJJvB4bAJ+L4Zx6MSxrCqPDavqcGwXc4Bb8Q607s44tRsbadjZQVtzD4lYCtt28Pk9hMI+CkuzKKrIoawyF59fxrEVJ65DJXdfIv1s3W+Avyul/jE4IQkhxPCjfGEUJoY/i/ycKE2bV2IYimgsgVIK19VYlkkw4MMwFFWj80gmbGo2NqNQTD65pN8J3vbVr/HANz9IV3M9F976Pea95xO71/V4TBKJFPFECtuOEwz69ro1vK/Sylwc26WxroPSytxD7jeVtNm+qYH6bS00bG+lub4d13XxeCwKSnOoGFvE+JNHkZUb6v+BO0qdbVEWL9zIile3smNry+7mfP6gl0DQi2kZJOIpot0JnMwQcKZlUDm2gJNOLuPkOVWUVeXJ8GfihHLQ5E5r/TPgZ0qpMaSTvIeAMqXUfwMPaq03DEqEQggxDASCfuyYD+WNkJ/dwtqaNQR8HpIpB43G7/Pi8ezdunXchHQfdDUbm9How9bgJWM9PP/Hb7L4vt+QXVLJzb9+mvJJs/cqo5TC7/fi8Vj0ROP09MTx+734fQeuqSqrSid0O2tbdyd32k1BtJFY0w5ad9YR62gn3hMnEY1Cog2jp4PiRDsFTpTOLpuNWxQrdnqpbcwhmvAxelI58y6YwjmXz2TM5LIBSZxaGrp44r7lLH2pBtdxqZ5QxCXXzGT81FKKy3MIRfbucsV1XFoau9m1o52tGxrZuLqOp+5fwZP/epOCkiymn1bNzDPGUFqZK4meGPH601q2Bvgu8F2l1FTSjSoeA47NeDZCCHEcCIV8OHE/ypdFbkSRSsRpr9tMYfXEg66jlNorwbNTLtOmlx3w9uiWZS/yyA9upa1+C6dedQsLbvkGvmDkoNs2TYNIOEA0liAeT6K1xu/z7Je4FJflYFoGdbWt6FPzcBtX4rTVYOBgJrqx6tfjadmGL96ISjSBzvTfp4BA+nVK0Z7t9agyVu8YxaP3rOdvP3+KyvHFvPOGM7jourmEsw/fgvdwUimHJ/65nIWPrkYZirMunsTpF0yguDznkOsZpkFhaRaFpVlMOzU9MkZXR4xVb2zjzVe38Oy/V/HMgyspqchhxhljmHnmGApLBn40kSOhtcZOOcRjqcwt5z1d56jM/zxeC3/Agy9wqEY64kSn+tNp5olg9uzZesmSJUMdhhBimIrHUqS2PkNs+V3Ur36BhxZ18q6v3sm086857Lpaa2q3tLL+rQaycgLMPHUUPn/6b+t4TyfP/u6rLHn4j+SVj+Xy//41Vaec2e+4tNbE4kmSSRu/z4Pf792vzO3fvJ+zZycYPypGtCvK6hdfI5xaT7bVgFJg+CJ48sekX9kVmJFizHAhpj8btIvT3YDbuplEwxoSLVtJtGwFrYlaVSx8azwPPuXH8gZYcPVsrr31AsqqC/odf1/129r46y9eoK62lTnnjecd184kJ//Y3ALu6ojx5qtbWfZyDVvWNQAwamwBs88ezYw55UQiFjhxtB3DSSZIxuLYyQR2MoXr2OC6aK0zL3BdAAXKSCfUhoEyTJQyUKaJYVoYpolhWZiWBYZJIu4S7bHp7k7S0ZagrSVGe2uMlqYorY1RuruTpJL9/51seUwCIS/hLD/hbD+R7ACR7ADhLD+RnACRbD/hrADBiI9wxI8/uH/yv38OoPde1rs8s57a3YXP3nq325ugppIOqZRDKmGTSjnYmddBKYVlGViWiWkZWB4z80rPszwmhjF8a1xdN32sVGZIwsGglFqqtZ59wGWS3KVJcieEOBTX1XRuWIhefy+dqx/lr890Me/aT3H+x77V72007upixbKdeDwG02aU07b5FR790SfpbKrjtGs+yTkfvK1f/dftS2tNNJYglXIIBnx4venEUSe7oWkFdvM64q07Wb9oEQXWJkyVwpNXTaBqDv6KUzELJqANL07KQcebIdqIEW/CiLdiuIn073LTB1YArCBusptY8xaiO1aRat+B9uaxqmUmf77PT2e3hwveM4frP3Uh5aMPPjTavsf2hf+8xaN/W0Ig5OP6j5/BlFlvf1xa7bpgxyDZAckusLvRqSjYMZxEDG3HUdoGO44ba8GNdeAmunCT0fQrlXm3E2gnhXZstGujHRtcG617a9T0gbv0zyRAuxOhvolRn3mqT1nUvuVVn/Kkk8i+5XfvX6dD0O7uRCz9e703MdOZ//pO68w0e5XTvdO983rXo0+Sd6CcYa9kRu3zUe0zeaDER+31ttf36/O2hz7ApN5nST/X3T27H7nQESRtvtKTyb/+vn6XfzsOldz1q587IYQ40RmGIql9hEJZGIaioLiA+g0rjmgbRSUR5p5RzZKX3uLer3yQ5rf+Q2H1RD74m2epmHxqv7ahtQYnAU4cnBRoG1yHAGDaDqmWHnTPVtzmdThtW0l0NNC+fSMeu4kSv4/AmDPxn3QxOnISyUSC7kQ7ntrleJP1eJLNvXvBNiJ0UkBS+7FMCNKDL9GG6q7DAEKhEKGpC0hEu+jetZmTk8/wkw/72BKbzR//tZCn73udy286ixs//w4iOQdPWNuau7nnVy+x6a16pp5aybUfPYNIduAwx8CFRCfEW9GJdnSyE1LdkIqh3AS4KXBt3GQPTlcTdncTdrQVJ9aJk+jCiXfjJHtwU8kDbl8ZJobHizI96do4w8QwTJQVQBnpWjgw9krM+gTX511nEsHeRMrNzHZ3l9G4mXwqvUz33hbvm4j1ft6ddOk9+929b7Vf8piZmyln7C6v9ks8Vaa4sSchy5TTysgs67PtPV92n4+ZuHoTLa33hKf3Tbv0PuvtdYIPkEjtM73P8t5vq/cqqthrcr8k8iDbPkDN5MEH5jrwfCc4tJ2GS3InhBD9FLc9ZPmzASgqiFCzedVeQ5D1x84VT7H8j5+hp72Z8nkfZMx5t6DCFbiu3u+2k9Ya7Cgku9CprkwtVBxwcVNxUm21pFpqSLXU4HTWY0dbcGIde5KHzG3DcE45gcJz8BWMxtA2RtPr0PQGvSmUY4bpcPLZ0VLAxh0e6jvC9CR9BENeAgEvqZRDV2cMpRRjRvs4fUqcYs92jJ4d+CyNv2IsqZLRdLfuYsz21/nO+2x2JKbw54da+MD9b/CBL7yTS288A7NPq16tNUte2sz9f3oN7Wqu+/iZzD1v/F7HUmsNiQ7oqYNYMzrRnk7i7Dh7kgONjreRaq8n2bkLu6sRO9qGE+vAScb3Op7K8mH4whi+XMzsMRi+fLS/ENdXgPbm4HhySRo5xJ0A0RjE4gbROMRiLilb46bASToYhsI0DQzLwO/3EAz7CIV9RLLTt0Nz8kOEI34M0xjWtxLFyHXY5E4pdQbwdaAqU14BWms9ZmBDE0KI4SXpeDD96YfwC3K8rFrVRFdLPVkFZYddt7O5jsd/9jnWvfQIJWNP5r0/fAB//njWrd7F6hV1bFzfSEVlLiWlEULeGDrRCvHWdC0d4KbiJJq3Eq9bRWLHUpyObfTW5hiWHzOYgxXMx587CsvrwfIGMX1BTG8QZXroTgbpTmXjyS6iOxmksd1iyw6HrTtdYql0R8tlFTmMmpzPrOp8SspzCIX3tEh1HJeWpm62b23h9ZpmttV4OHXWdGaM6cTXswmrawe5RaVEci+gp72ZUTtXc9u1q9nZWcZ9f9jMI395iU9+5xpOOWM8rU3d/OsPr7B2+Q5GTyjifZ84m4KSrPTzbPFWaN8M3TvTyVxvTRYKTA+uNki0NZBsXEeqdQt2VwNOomd3nMryo335OOFJqMgoAqWTCFacTGDU9PQzhEKcAPpTc/cn4DPAUuAQT0MKIcTIZms/hjcEhkl+VrqlYsPGVYdM7lzXZdkjd/LM776CYydZcMs3OO3a/8K00l2XzDtrNM1NPezc2gBdtShvEtfr4mpFtKOTri1LSW17CbNzLUq7KNPCEykiWDEFTyCCJ5iD4Q2gDS9xcui2IySChYSLy3llcSvrNiXoSvhwnHQNo8rUJHk9JuVVRZx6Tj6V1fmUjsrF49lTs+a4LrFYkmg0SSJuk0w6pBIOwVCAiVMrmDi1gq7OOM+vDhPx5zBt/Gyy3Fqs9hqyfSHCeUXEmrdjNW7hU5cvprFzHY9/bzn3Zp9PwvFjWgZX3TyXsy6cgIo3oneuhK5taDuWDsD0ob05pFyTeHMNybqVOE1vYXfV9bZmwPBnk/JXYuePJ3fSWWRPPBszMjBdswhxPOlPctehtX58wCMRQohhTlvBdILkDZObaZRav+FNxp920QHLN9Wu59EffpJtq16hevrZXPqFX5JfMXbvQnaMfE8d+RVNaKCnx2XzspcwdzyBJ7oFpTSBYC6+son4sgvxhvLQngjdqphmOx/bV0LxhCkE84qJGAYR0rczu7pjzH6nZlJXgs72GM88tJL2lh6u//hZ5OWHyM4N7teVhuO6pFI2yaSN42QSKBMCIYtAaM+vi2TCobsrQTJlkZUXJtrtZfE6TU54EuMqxpLracXsrCUcKiJUNIZoyzas+o28N+s1Es4qNrRNxJNdQWViB6mlQUyvB1C4GCSTDsnuZlKN63Ca1uH0NKV3qhRmqAAnbxZW2TxKTr8WK7f66E+qECNQf5K755VSPwQeABK9M7XWywYsKiGEGIYMjxfXNlFWCIskpWOnsmbhg5x14xf3qi2K93Ty8t0/5tV//hJvIMTl//1bpl9yw97PkzkpdFctRBtAKdqb26h/+V7CHYsIaBtPKJdA5VT8ueUYoWK6rErq7HLa3DK6E2GSicz4sj2w9Y1O8gsdyiqyKSgKYxiKUNBPV3eM7LwgZRW5NO3o4IE7XiMnO0BeQXhPHDo9hFoiaWNnxsrt6U7S0thDMmbTuKONHVuacVIu/qCH/OIIheXZlFXmUjk6B0j3TbdrezsNOzp4emEKAx9TJp/M2PLJeJ1WQrnbCZZOIdlcQ9eu9ZziWY52l8I6aFy/z0HWe461GSokWXQ6unAOFWddg5UzSmrlhOiH/iR3czPvfZvbamD+sQ9HCCGGL6/XwnH8aE8EJ17H7Iuu45HffIvX7/8tp77rY3Q17WTFk3/n9ft/S097EydfcB3nf/zbRPJLdm9Daw3RXenEznVoqW+k4YXbyUqsItu0CBZVESysRuWMpVmNZ0N3OYWlJ1E5roi8PuPIOrZLLJaioz1GW2uUpl3dNNR34g94GDOugIqqHAJ+L7F4kng8xYRp6VvH61fWUVCShetqkimbZDKF62q0q2ms76ZuRwcGipWvbGbzW/WYlkFpZS45eSG8Xot4T4q6za3s2tpGIOwltyBMVl6Q/KIIo8YW7G7EqF1NXENDczOOqiCv1CZc1kL+6Fp0zy7srqZ069VUHNex0coipUJEdQF27hSKTp5PpKRKkjkh3ob+jFBx3mAEIoQQw53fb+HGvShvBLcrzslnnMGaN87niV9+kSd/9T+7+z8bM2s+82/5GuUTZ+21vnYS6PZNkGijp7OHLU/9gdye18j1eAlVTCZYOpXu7Nms6RlPceVoyirzKT9IcmNaBuGIj3DER/moHNxpmubGbmo2NrNmVT07trUx5ZRSvF6LRDJFMMfH6RdNIB5P0hONk8p0KKs11G3rYOe2drw+ix0bGlj6wiZyC8O899azmH7aaLy+/nWsoLXGdlwc20m/Ow45JUVAEWhoS6SwfafiG+XFMEwSCZuWxm60qykqyaawKLL7mUAhxNt30CtWKXWD1vpupdRnD7Rca/2TgQtLCCGGn0DQi9Phw/BHsO0UKtrM9f/3L1Y/ex8t2zYQzClg3NwLKagcv9+6OtGJbluH6yTZtOhpglv/Sh4OobIJ+KtOpzXnfOpUNeMmlzPTc+S9VBmGoqgkQmFxmF11naxb3cCrL22lqjqXUaNzsR2H8999CmhNKumAUuzY0s7O7e14vSaxzhhPPLIKj8fk0vfN5px3TMbjPbI4lFJ4LBNPpssTrTWuq7EdB9fVeLwm6HQXIZZlkpNjUFKSc8TfVQhxaIe6cnvHfDn44IZCCHECCQS8xJQfbzAbG3C6duC1PJxy0fWHXE9HG9Dtm0nEoux45NuEYpvwZRcRHr+A9rJ30+4fy6jqfMqPwVihSilKy7PJLwyzcW0jtVta2VbbRlFxBK/XYPHCjeSXZOMLpLs/UY7D0/9YTqwnyWkLTuLia2aSlXPoToSPJBbTVDIGqhCD7KDJndb695n3bwxeOEIIMXxZloGjfPgjEaKA07XrsOvo7p3ozi20ba+h56XvEXB6yBozh+TED9GQNZey8hx8R1hD1h9er8mUU0qpHptH7ZY2Gnd1kYinKKkqoKstys7NzWxdW08ilmLyzFFcdsNsSkflHvM4hBCDrz+dGI8Bfg7MI92Q4lXgM1rrmgGOTQghhhXDUDh4Mf3pGxpud+Mhy+uuHeiurTStfYPUkp9jebxEpl1L67hbycovojo7MOANBkJhH5OnlTB5Wsnu7k3efHULdk+csy6exCnzqhk1pmBAYxBCDK7+/Ln4N+DXwFWZ6euAv7OnFa0QQpwQlFIkHQvDl+5KxIm1HbSsjjakE7sVz5BceQeWP0Lg1P+ipexdFBVnEwx4Byvs3Xpvj846cyyzzhx7mNJCiONVfx6ECGqt/6q1tjOvuwH/QAcmhBDDUdL1YHgDYBg48c4DltGJdnT7JtrXvURq5R14grn4z/wqnaPeQ1lZ7pAkdkKIE0d/au4eV0r9D3Av6duy1wKPKaXyALTWrQMYnxBCDCu29mIoA+UJ48a7ce0EhrVnDFbtJNBt64huW05s6e2Y/gi+M75Md+H5lJZk7TXElxBCDIT+JHfXZN4/us/860gne2OOaURCCDGMOcoHWqGtEE4qlh5hIqsSSHf9odvW43Rso2vxb0EZBOZ+lo6CBZQWS2InhBgc/enEePRgBHIsZYZLuwxIApuBm7XW7UMalBBiRFCGF+0aOCqMk2yA6K7dyR3d26G7nvaXfoMd7SA86ybaSy6luCiC1yuJnRBicBz0mTul1KlKqZI+0zcqpf6tlPpF7y3ZYexpYKrW+mRgA/ClIY5HCDFCmJaJo3ykjAhuKg49DQDoRAe6s5auVf8m3rSJUPVcesZ9mJzcLAJ+zxBHLYQ4kRyqQcXvSdd8oZQ6G/g/4C9AB3D7wIf29mmtn9JaZ0bV5jWgYijjEUKMHF5fOrlTvixcO4nu2pl5zm49dtN6opsXYgay0NM/iRnMIztL2p8JIQbXoZI7s09jiWuB27XW92utvwKMG/jQjpkPAo8faIFS6hal1BKl1JKmpqZBDksIcTzyeS0c5cMXyQfAbq1BNy6HWDM9q+7HjnUSmno1XcGpFOSHZOB7IcSgO2Ryp5TqfSZvAfBcn2XHvjv1I6SUekYptfoAryv6lLkNsIF7DrQNrfXtWuvZWuvZhYWFgxW6EOI4Fgh6cJSPQF4xKEWqfTtol551z9JdtwZ/8UnERl9HdnZIGlAIIYbEoZK0vwMvKKWagRjwEoBSahzpW7NDSmt9/qGWK6VuAi4FFmit9aAEJYQY8fx+Dz3KR1Z+PjEM7O4m3C3PENv6Yrp17Cnvp9VbQpHcjhVCDJFDjS37HaXUs0Ap8FSfBMkAPjkYwb1dSqmLgS8C52ito0MdjxBi5DBNE0f5MS2LZHAc8Y4mzJ1vEm+vJ3vShXTknk12dhDDkNuxQoihccjbq1rr1w4wb8PAhXPM/ArwAU9nnnd5TWv9saENSQgxEpimImWEQUG7NQ1v5wO0bnodb24F6qSrcbx5ZEV8h9+QEEIMkCF/dm4gaK2PpwYfQojjiFKKhA4CsLVjLG2tE6gssymeeQVNgalEwn5pRCGEGFL9GVtWCCFEH0nHh6tNLN3D42vPYMyCG+gOTcMxw1JrJ4QYciOy5k4IIQaSYRh06XwueYfD2QsseuIuPcVTCQW9mKb8zSyEGFryU0gIIY5QIOih0zMGwzTwWQ7b7am4eIhEpIWsEGLoSc2dEEIcoVDYT2uqhF/fk49ja67+0Lh0oifjxwohhgFJ7oQQ4gj5vCZKQTgni1nzRuO4muxsnzSkEEIMC5LcCSHEEfJ6LZRSXH3jXGzbpas7QTjkHeqwhBACkOROCCGOmGkqLMugozOO62ppSCGEGFbkp5EQQhwhpRSRsA/H0SilyMkJDHVIQgixm9TcCSHE25Cd5cfntfB4DCxLGlIIIYYPSe6EEOJtUEoRCHiGOgwhhNiP3JYVQgghhBhBJLkTQgghhBhBJLkTQgghhBhBlNZ6qGMYFpRSTUDtIOyqAGgehP2I/pNzMjzJeRl+5JwMT3Jehp/BOCdVWuvCAy2Q5G6QKaWWaK1nD3UcYg85J8OTnJfhR87J8CTnZfgZ6nMit2WFEEIIIUYQSe6EEEIIIUYQSe4G3+1DHYDYj5yT4UnOy/Aj52R4kvMy/AzpOZFn7oQQQgghRhCpuRNCCCGEGEEkuRskSqmLlVLrlVKblFL/M9TxnCiUUqOUUs8rpdYopd5SSn0qMz9PKfW0Umpj5j03M18ppX6ROU8rlVIzh/YbjGxKKVMptVwp9WhmerRSanHm+P9DKeXNzPdlpjdlllcPaeAjmFIqRyl1n1JqnVJqrVLqNLlehpZS6jOZn1+rlVJ/V0r55VoZfEqpO5RSjUqp1X3mHfG1oZT6QKb8RqXUBwYiVknuBoFSygR+DVwCTAauV0pNHtqoThg28Dmt9WRgHnBr5tj/D/Cs1no88GxmGtLnaHzmdQvw28EP+YTyKWBtn+nvAz/VWo8D2oAPZeZ/CGjLzP9pppwYGD8HntBaTwROIX1+5HoZIkqpcuC/gNla66mACVyHXCtD4c/AxfvMO6JrQymVB3wNmAvMAb7WmxAeS5LcDY45wCatdY3WOgncC1wxxDGdELTW9VrrZZnPXaR/UZWTPv53ZYrdBVyZ+XwF8Bed9hqQo5QqHdyoTwxKqQrgncAfM9MKmA/clymy73npPV/3AQsy5cUxpJTKBs4G/gSgtU5qrduR62WoWUBAKWUBQaAeuVYGndb6RaB1n9lHem1cBDyttW7VWrcBT7N/wnjUJLkbHOXA9j7TOzLzxCDK3J6YASwGirXW9ZlFu4DizGc5V4PnZ8AXATcznQ+0a63tzHTfY7/7vGSWd2TKi2NrNNAE3Jm5Xf5HpVQIuV6GjNZ6J/AjYBvppK4DWIpcK8PFkV4bg3LNSHInTghKqTBwP/BprXVn32U63WRcmo0PIqXUpUCj1nrpUMci9mIBM4Hfaq1nAD3suc0EyPUy2DK37K4gnXiXASEGoKZHHL3hdG1Icjc4dgKj+kxXZOaJQaCU8pBO7O7RWj+Qmd3Qe/so896YmS/nanCcAVyulNpK+jGF+aSf9crJ3HqCvY/97vOSWZ4NtAxmwCeIHcAOrfXizPR9pJM9uV6GzvnAFq11k9Y6BTxA+vqRa2V4ONJrY1CuGUnuBscbwPhM6yYv6YdhHx7imE4ImWdN/gSs1Vr/pM+ih4HeVkofAP7dZ/6NmZZO84COPlXu4hjRWn9Ja12hta4mfT08p7V+H/A8cHWm2L7npfd8XZ0pPyz+Qh5JtNa7gO1KqQmZWQuANcj1MpS2AfOUUsHMz7PecyLXyvBwpNfGk8CFSqncTK3shZl5x5R0YjxIlFLvIP2MkQncobX+ztBGdGJQSp0JvASsYs+zXf9L+rm7fwKVQC1wjda6NfPD81ekb3tEgZu11ksGPfATiFLqXODzWutLlVJjSNfk5QHLgRu01gmllB/4K+lnJluB67TWNUMU8oimlJpOupGLF6gBbiZdESDXyxBRSn0DuJZ06//lwIdJP6cl18ogUkr9HTgXKAAaSLd6fYgjvDaUUh8k/XsI4Dta6zuPeayS3AkhhBBCjBxyW1YIIYQQYgSR5E4IIYQQYgSR5E4IIYQQYgSR5E4IIYQQYgSR5E4IIYYJpVS+UurNzGuXUmpn5nO3Uuo3A7C/K2WcayFGHmktK4QQw5BS6utAt9b6RwO4jz8Dj2qt7ztcWSHE8UNq7oQQYphTSp2rlHo08/nrSqm7lFIvKaVqlVLvUkr9QCm1Sin1RGZEFpRSs5RSLyilliqlnuztRb/PNk8HLgd+mKkdHDv430wIMRAkuRNCiOPPWNJDtl0O3A08r7WeBsSAd2YSvF8CV2utZwF3AHt1nK61foV0L/pf0FpP11pvHswvIIQYONbhiwghhBhmHtdap5RSq0iPevNEZv4qoBqYAEwFnk53lI8JyLBgQpwgJLkTQojjTwJAa+0qpVJ9xg51Sf9cV8BbWuvThipAIcTQkduyQggx8qwHCpVSpwEopTxKqSkHKNcFRAY1MiHEgJPkTgghRhitdRK4Gvi+UmoF8CZw+gGK3gt8QSm1XBpUCDFySFcoQgghhBAjiNTcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIJLcCSGEEEKMIP8f/4JPmekRD3sAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test.result.plot_spin_trajectories(plot_type=\"spins\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAACZCAYAAABNGxy5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaLElEQVR4nO3de5hddX3v8fdn77nlwi0khEtIQjAqF4ViygHtoZzCU9BSqJZWqlRRFDhHT7H2gFDaU2yxPhaK5ZRqpQjaIl6eoGBRFII+SltRE4KEawkgEJpASIAkhJnZM/t7/lhrTybJXPbM7LX2nrU/r+eZJ3td9p5vZs3a+zu/3+/7+ykiMDMzM7NiKDU7ADMzMzNrHCd3ZmZmZgXi5M7MzMysQJzcmZmZmRWIkzszMzOzAnFyZ2ZmZlYgHc0OoFXMnTs3Fi9e3OwwzMzMzMa1atWqFyNi3kjHpm1yJ+n3gMuBw4BjI2Jlun8x8AjwWHrqvRFxwXivt3jxYlauXJlNsGZmZmYNJOnp0Y5N2+QOeBB4F/CFEY49ERFH5xuOmZmZWfNN2+QuIh4BkNTsUMzMzMxaRlELKg6RtFrSjyT992YHU1RPrXmaS069gotOupwf3HxPs8MxMzMzWrzlTtIKYP8RDl0WEbeN8rT1wMKI2CTpLcCtko6IiC0jvP55wHkACxcubFTYbWPd4xvYsmkrHZ1lnnrw2WaHY2ZmZrR4chcRJ0/iOX1AX/p4laQngNcDu1VLRMR1wHUAy5Yti6lF234qvf0AzN5nNtWBwSZHY2ZmZlDAbllJ8ySV08dLgKXAk82Nqpj6eysAzJjdw6CTOzMzs5YwbZM7Se+UtA44HviOpO+nh04AHpB0P7AcuCAiNjcpzEKr9FWQRFdPJ4MD1WaHY2ZmZrR4t+xYIuJbwLdG2H8LcEv+EbWfSt8And0dlMtlBged3JmZmbWCadtyZ83X39tPV08X5Y6Sx9yZmZm1CCd3Nmm1lrtSueQxd2ZmZi3CyZ1NWn9fhc7uTsodZQar7pY1MzNrBU7ubNIqvf109STJXdUFFWZmZi3ByZ1NWn9vha7uLnfLmpmZtRAndzZplb5KUi3bUXK1rJmZWYtwcmeTVukfoLO7k1LZ1bJmZmatwsmdTVp/X2VozJ1b7szMzFrDtJ3EeLp5/umNfPFPv9LsMBpq47ObWHTYgnTMnZM7y9fNf/1NfvnQMwDste+enP+376Oj029pZmZ+J8xJuaPEnPl7NzuMhpqz/z4c+45jWH33GhdUWO5Wfv9+Zu8zixmze3joJ4/x8sYtzD1wTrPDMjNrOid3OZl70L5ccPU5zQ4jEw/86GGPubPcVatVjj7xSBYevoAbLruZSl+l2SGZmbUEj7mzKSt3uFvW8hfVQCXR1dMJJFPzmJmZkztrABdUWDNUB6uUyiU6u5Pkzi13ZmYJJ3c2ZZ4KxZqhWg1KJdHZnYwuqfQNNDkiM7PW4OTOpqzWchcRzQ7F2kREEBFuuTMzG8G4yZ2k/y1pnzyCsemp3JH8GlWr7pq1fNR+11Qq0dXtMXdmZsPV03I3H/i5pG9IOlWSsg7KppdSOfk1clGF5aWajvEslUTnUEFFfzNDMjNrGeMmdxHxZ8BS4IvAOcDjkv5a0qEZx2bTRLmjDOC57iw3UU2GAJTKpaFqWY+5MzNL1DXmLpLBVBvSrwFgH2C5pL/JMDabJmotd1VXzFpOqmlyp5I85s7MbBfjTmIs6ULgfcCLwPXARRFRkVQCHgcuzjZEa3VuubO81f6QKJd3jLlzcmdmlqhnhYo5wLsi4unhOyOiKum0bMKy6cRj7ixvteROpRLljjKlktwta2aWqie5uwZA0vBFG7dGRCUiHskmLJtOatWybrmzvNSqZUulpL6rq6fLBRVmZql6krv7gIOBlwABewMbJD0PfDgiVmUXnk0HtW7Zm6+4ha4ZXcyY3cMfXPpOunq6mhyZFdVQtWzaatzZ08kvfvQQLz63eafz/LtoZu2onoKKu4B3RMTciNgXeDtwO/C/gM9lGZxND4uPOJjFRx7Mtlde5bm16/npd+9j/ZPPNzssK7Dh1bIAv3rK0cyY3cPmDS8Nffl30czaVT0td8dFxIdrGxFxp6SrIuJ8Sd0ZxmbTxPxF87j4xo8C8PBPHuPaP7rB4+8sU4NDY+6SbtkzP/7bu53z4L8/yuc+duNQZa2ZWbuoJ7lbL+kTwNfS7XcDz0sqA/4Et524ctbyUGu5K5dH73yojccLJ3dm1mbq6ZZ9D7AAuBX4Fsn4u/cAZeD3M4tsHJKulPSopAckfUvS3sOOXSppraTHJJ3SrBjbkZM7y0MMW35sNJ5/0cza1Zgtd2nr3DUR8d5RTlnb+JDqdhdwaUQMSPoMcCnwCUmHA2cBRwAHAiskvT4inG3kwB+olofBXQoqRlJL/Ab9u2hmbWbMlrs0IVokqeVKzSLizoioTWx1L0nrIsAZwNcioi8iniJJQI9tRoztaMe0KP5AtewMFVSURl/qeke3rH8Xzay91DPm7kng3yV9G3i1tjMirs4sqon7IPD19PFBJMlezbp0n+XA3bKWh+ouBRUjqf0uuhXZzNpNPcndE+lXCdgj23B2JmkFsP8Ihy6LiNvScy4jWe/2K5N4/fOA8wAWLlw4hUitxh+olofqLlOhjKSW+Lla1szazbjJXUR8EkDSzIjYnn1IO33vk8c6Lukc4DTgpIiovYM/R1L0UbMg3TfS618HXAewbNkyfwI0wI6lyNxyZ9mJoRUqxq+W9R8aZtZuxq2WlXS8pIeBR9PtoyQ1ffJiSacCFwOn75J0fhs4S1K3pEOApcDPmhFjO6qNufMHqmVp1xUqRlI75qlQzKzd1NMt+3fAKSRJExHxC0knZBlUna4FuoG7JAHcGxEXRMRDkr4BPEzSXfsRV8rmZ8eYOyd3lp1dJzEeSa1atuqCCjNrM/Ukd0TEs2kCVdP0ZCkiXjfGsU8Bn8oxHEu5W9byUE+1bNnT8phZm6onuXtW0luBkNQJXAg8km1YNl25oMLyUE+3rAsqzKxd1ZPcXQBcQzKdyHPAncBHsgzKpq+itdy98MxGfrz83kknq0f9jyN5w7JDGxyVDVXLuqDCzGw39VTLvgiMtkKF2U6KNonxvbffxw+++m/M2nPGhJ/72rZeXnjmRSd3GXBBhZnZ6MZN7iTNAz4MLB5+fkR8MLuwbLoq2iTG/b399Mzs5sq7L5/wc//23M+51SgjQ1Oh1LH8mAsqzKzd1NMtextwD7CCFiiksNZWtG7Z/t4KnT2dk3puqVxycpeRWrfsWNWyXufYzNpVPcndzIj4ROaRWCGUSiVKJRXmA7XSV6Gre3LJnUpO7rJS+7mWx+qWHVpb1t2yZtZexp3EGLhd0jsyj8QKo1QuFWbMXaWvQmd3XTMG7aZISW6rqaflTi6oMLM2VU9ydyFJgveapC2StkraknVgNn2VO8qF6Zat9A3QOcmWu1K5RDWcWGRhqKBijGrZoWl53HJnZm2mnmrZPfIIxIqjSGPN+nv7p9AtK3cJZqSegood3bLF+F00M6vXqO+Mks4e9vhtuxz7aJZB2fTmlrtEuVwqzM+h1dTXLZsW9xTkDw0zs3qN1S378WGP/36XY54GxUZV7igV5gO1v6+frklWy6pUcstdRjzPnZnZ6MZK7jTK45G2zYaUyiWqhSmomMKYu5I83isjO8bcjTEVigsqzKxNjZXcxSiPR9o2G1Kkbtn+3n66erom9dwijT1sNTu6ZcdvufM1MLN2M1ZBxRslPUDSSndo+ph0e0nmkdm0VSoXp1s2abmb3FQoSUFFMX4Oraaeee6ktOXOradm1mbG+tQ6LLcorFDKHWWqBWm5S+a5m/xUKEVJcltNLWkes6BC8lyDZtaWRk3uIuLpPAOx4ih3lNjw9EZW3PTjZocyZVNK7lxQMapqtcq9t69i+5bXJvX8taufAsYuqIB0lRC3nppZm5lcf5PZGOYtmMvqH6zhm9d8p9mhNMT8RXMn9Ty5oGJU6598npv+avmUXmOPObPpnjH2eMiS5xo0szbk5M4a7txPv4f+3kqzw2gICbpndE/quS6oGN1A/wAAH/r0ezn8rW+Y1Gt0dnUMrUIxGl8DM2tHTu6s4UqlEj0zJ5cQFUm5XHJBxShqaw/3zOrO9HfFradm1o7GTe4krWH3qU9eAVYCV0TEpiwCM5vuVHKr0WhqU+WM1/I2VWW33JlZG6qn5e4OYBC4Od0+C5gJbAC+BPx2JpGZTXOlcsmtRqOoJXfjFURMVbJKiJM7M2sv9SR3J0fEMcO210i6LyKOGb7+rJntzIP5R1frls265c4Jtpm1o3r+bC5LOra2IelXgdo78kAmUZkVgOe5G93QJMQd2bbcuaDCzNpRPS13HwJukDSbZHWKLcCHJM0CPp1lcGbTmTyB7qjyGnPn1lMza0fjJncR8XPgTZL2SrdfGXb4G1kFZjbdlTzea1Q7Wu6yTe5UEoODxVgtxcysXvVUy3YDvwssBjpq6zVGxF9mGpnZNFcb7xURQ+ucWiKvgopkOhq33JlZe6mnW/Y2kqlPVgF92YZjVhyl0tAfQk7udrGjoCL7alkXVJhZu6knuVsQEadmHskESbqSZBqWfuAJ4AMR8bKkxcAjwGPpqfdGxAXNidLaWa1VqjpYpVTKNomZbmrdslm33LmgwszaUT3vrP8h6U2ZRzJxdwFHRsSbgf8ELh127ImIODr9cmJnTaG05c4tR7vLt6DCyZ2ZtZd6Wu5+DThH0lMk3bICIk2qmiYi7hy2eS9wZrNiMRtJrbXOLUe7y7Og4qUXXuG+FQ/Udf4BS+ZzwJL5rL3/Kba8uHXc8/feby+WvHnRVMM0M2uoepK7t2cexdR9EPj6sO1DJK0mmbblzyLinpGeJOk84DyAhQsXZh6ktZdal6NbjnaXV0HF7L1n8+jPHuf6S79S1/n7LZzLRTd+hM+e9wUixm9xLZXEVT/8pNdSNrOWMmpyJ2nPiNgCjP/na0YkrQD2H+HQZRFxW3rOZSSTKdfevdcDCyNik6S3ALdKOiL9v+wkIq4DrgNYtmyZ+86sodwtO7q8CirOv+oP2bz+5brO/e71d/OfK9fy2rZeIoLT/+cpHHXiEaOef9+KB/jOP62g99VeJ3dm1lLGarm7GTiNpEo2SLpjawJYkmFcyTeJOHms45LOIYnxpEj/zI6IPtKq3ohYJekJ4PXAymyjNdtZuexu2dHkNeaue0Y3ByyZX9e5e83dg0rfAJW+CgBzF+w75nPnHjQHgEqfF+oxs9YyanIXEael/x6SXzj1k3QqcDHw6xGxfdj+ecDmiBiUtARYCjzZpDCtjdVa7jzP2u7y6padiK6eLvp7++nvraTbnWOe39mdHO/v7c88NjOziRj3nVXS29KlxpB0tqSrJbXCALVrgT2AuyTdL+kf0/0nAA9Iuh9YDlwQEZubFKO1MaUFFV5fdnfJ9DBqqSliuno6qVaD3leT6TxrydtY54Nb7sys9dRTUPF54ChJRwF/AlwP/Avw61kGNp6IeN0o+28Bbsk5HLPdlF1QMarBgWpLtdrBjmRu+5bt6fbYb4+182vduGZmraKed9eBdDzbGcC1EfEPJC1mZjaGoYIKt9ztpjpYzXy83UTVkrltL9eSu3q7ZZ3cmVlrqaflbqukS4GzgRMklYCx3/XMbKcVKmxngwODLZfcdfV0AfDqK9vT7fq6ZT3mzsxaTT0td+8mqT49NyI2AAuAKzONyqwAhpI7F1TsZnBgsAW7ZZO/dXd0y9bXcucxd2bWasZtuUsTuquHbT8D/HOWQZkVQcnVsqMaHKxmPsfdRNWStYm23HnMnZm1mrEmMd5KMp/dbodIlh/bM7OozAqg1jLlatndVVuwoKKWrNWSu/pb7pzcmVlrGWueOxdNmE1BbZoPV8vubnBwGoy5Gze5S94+XVBhZq2mtf50NisQV8uObnCgdatlt299jXK5NG58brkzs1ZVT7WsmU1CrdvxubUbqGMN+raydfO2luuWrSVrWzZtHbdLFpKW2Y7OMhvXbeKXDz2bdXgTNnOPHvZbOG9Sz630V/ivBv/e9szqZv/F+zXuBc1sVE7uzDLSMytZTP6mv1re5Eha05I3LWp2CDuZuecMIOmW3ffAfep6zqy9ZvKzO1bzsztWZxnapF3xr5cwZ//6/i/D/evn72TFTT9ueDx//vWP173Wr5lNnpM7s4wsefMi/vgL59P3mudBG8mBh7bWh/yec/bgE1/+KFtfepX9Dt63rud87B/PZ+O6TRlHNnG/fPBZvnv9Cra+9Oqkkrstm7ayx5zZ/OH//b2GxPNfazdw67V3sPWlbRxAa113syJycmeWkVKpxNJjljQ7DJuARYcfPKHz5y+ax/xFk+v6zFKty3uy4wErfQPM2msmR77tjQ2JZ1baKuo5Ac3y0VqDXszMbMp2zME3uWSqv7d/3Hn+JmLHUm1uxTbLg5M7M7OCmWoyVekboLOrccndjqXaXFlslgcnd2ZmBdM1xaXRKv2VuiqG69WZziHoaWPM8uHkzsysYDp7ptpyV6Gru3FDsmuv5TF3ZvlwcmdmVjBTH3NXGWptawRP+GyWLyd3ZmYFM9Vkqr+30tCCio6uWsudkzuzPDi5MzMrmK4pJneVvsaOuSuVSnR2dbigwiwnTu7MzAqm3FGmVNKkk6lkzF3jkjtIuoqd3Jnlw8mdmVkBdfV0TarlLiLSMXeNTe46uzvdLWuWEyd3ZmYF1DnJlrKBSlKE0chu2drrObkzy4eXHzMzK6Cu7k62b9nOlk1bJ/S817b1ps9v7MdDZ3cH27e8NuF4zKajzu4OZsye0bTv7+TOzKyAemb3cN/da7jv7jWTfn4jzZjVw0M/eYxLTr2ioa9r1opO+N3jOOuSdzbt+zu5MzMroLP//EyeeXjdpJ5b7ixzzMlvbmg8777knTz1wNMNfU2zVnXAkvlN/f5O7szMCmjRYQtYdNiCZocxZMHSA1iw9IBmh2HWFlxQYWZmZlYgTu7MzMzMCsTJnZmZmVmBKCKaHUNLkLQRyGO071zgxRy+j9XP16Q1+bq0Hl+T1uTr0nryuCaLImLeSAec3OVM0sqIWNbsOGwHX5PW5OvSenxNWpOvS+tp9jVxt6yZmZlZgTi5MzMzMysQJ3f5u67ZAdhufE1ak69L6/E1aU2+Lq2nqdfEY+7MzMzMCsQtd2ZmZmYF4uQuJ5JOlfSYpLWSLml2PO1C0sGSfijpYUkPSbow3T9H0l2SHk//3SfdL0n/L71OD0g6prn/g2KTVJa0WtLt6fYhkn6a/vy/Lqkr3d+dbq9Njy9uauAFJmlvScslPSrpEUnH+35pLkl/nL5/PSjpq5J6fK/kT9INkl6Q9OCwfRO+NyS9Pz3/cUnvzyJWJ3c5kFQG/gF4O3A48AeSDm9uVG1jAPiTiDgcOA74SPqzvwS4OyKWAnen25Bco6Xp13nA5/MPua1cCDwybPszwGcj4nXAS8C56f5zgZfS/Z9Nz7NsXAN8LyLeCBxFcn18vzSJpIOAPwKWRcSRQBk4C98rzfAl4NRd9k3o3pA0B/gL4L8BxwJ/UUsIG8nJXT6OBdZGxJMR0Q98DTijyTG1hYhYHxH3pY+3knxQHUTy8/9yetqXgd9JH58B/HMk7gX2luTVzjMgaQHwW8D16baA3wCWp6fsel1q12s5cFJ6vjWQpL2AE4AvAkREf0S8jO+XZusAZkjqAGYC6/G9kruI+DGweZfdE703TgHuiojNEfEScBe7J4xT5uQuHwcBzw7bXpfusxyl3RO/AvwUmB8R69NDG4D56WNfq/z8HXAxUE239wVejoiBdHv4z37ouqTHX0nPt8Y6BNgI3Jh2l18vaRa+X5omIp4DrgKeIUnqXgFW4XulVUz03sjlnnFyZ21B0mzgFuBjEbFl+LFISsZdNp4jSacBL0TEqmbHYjvpAI4BPh8RvwK8yo5uJsD3S97SLrszSBLvA4FZZNDSY1PXSveGk7t8PAccPGx7QbrPciCpkySx+0pEfDPd/Xyt+yj994V0v69VPt4GnC7plyTDFH6DZKzX3mnXE+z8sx+6LunxvYBNeQbcJtYB6yLip+n2cpJkz/dL85wMPBURGyOiAnyT5P7xvdIaJnpv5HLPOLnLx8+BpWl1UxfJYNhvNzmmtpCONfki8EhEXD3s0LeBWpXS+4Hbhu1/X1rpdBzwyrAmd2uQiLg0IhZExGKS++EHEfFe4IfAmelpu16X2vU6Mz2/Jf5CLpKI2AA8K+kN6a6TgIfx/dJMzwDHSZqZvp/VronvldYw0Xvj+8BvStonbZX9zXRfQ3kS45xIegfJGKMycENEfKq5EbUHSb8G3AOsYcfYrj8lGXf3DWAh8DTw+xGxOX3zvJak22M78IGIWJl74G1E0onA/4mI0yQtIWnJmwOsBs6OiD5JPcC/kIyZ3AycFRFPNinkQpN0NEmRSxfwJPABkoYA3y9NIumTwLtJqv9XAx8iGafleyVHkr4KnAjMBZ4nqXq9lQneG5I+SPI5BPCpiLix4bE6uTMzMzMrDnfLmpmZmRWIkzszMzOzAnFyZ2ZmZlYgTu7MzMzMCsTJnZlZi5C0r6T7068Nkp5LH2+T9LkMvt/veJ1rs+JxtayZWQuSdDmwLSKuyvB7fAm4PSKWj3eumU0fbrkzM2txkk6UdHv6+HJJX5Z0j6SnJb1L0t9IWiPpe+mKLEh6i6QfSVol6fu1WfSHveZbgdOBK9PWwUPz/5+ZWRac3JmZTT+HkizZdjpwE/DDiHgT8BrwW2mC9/fAmRHxFuAGYKeJ0yPiP0hm0b8oIo6OiCfy/A+YWXY6xj/FzMxazB0RUZG0hmTVm++l+9cAi4E3AEcCdyUT5VMGvCyYWZtwcmdmNv30AUREVVJl2NqhVZL3dQEPRcTxzQrQzJrH3bJmZsXzGDBP0vEAkjolHTHCeVuBPXKNzMwy5+TOzKxgIqIfOBP4jKRfAPcDbx3h1K8BF0la7YIKs+LwVChmZmZmBeKWOzMzM7MCcXJnZmZmViBO7szMzMwKxMmdmZmZWYE4uTMzMzMrECd3ZmZmZgXi5M7MzMysQJzcmZmZmRXI/wfXmGw4Pg/38AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test.result.plot_spin_trajectories(plot_type=\"energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GPU Acceleration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now showcase our implementation with GPU acceleration using Pytorch and CUDA. While the speedup advantage is minor for smaller problems (such as an $N=20$ instance), the difference is more apparent at large spin counts. This is done by specifying the argument use_GPU=True, and requires a CUDA installation. By default, the argument use_GPU is set to False, and will run the solver in batches in parallel on the CPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No External Field Detected\n",
      "Target Ising Energy: -inf.\n",
      "Best Ising Energy Found: -27.0.\n",
      "Corresponding Spin Configuration: [-1. ...  1.].\n",
      "Time Elapsed: 0.7007176876068115.\n",
      "Number of Runs Completed: 1.\n"
     ]
    }
   ],
   "source": [
    "test = Ising(J).solve(cac_gamma=gamma, use_GPU=True, hyperparameters_randomtune = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To specify the number of solver runs to be executed in parallel on a given device, the argument num_parallel_runs (set to 1 by default) can be modified. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No External Field Detected\n",
      "Target Ising Energy: -inf.\n",
      "Best Ising Energy Found: -29.0.\n",
      "Corresponding Spin Configuration: [-1. ... -1.].\n",
      "Time Elapsed: 0.7409710884094238.\n",
      "Number of Runs Completed: 5.\n"
     ]
    }
   ],
   "source": [
    "test = Ising(J).solve(cac_gamma=gamma, use_GPU=True, num_runs=5, num_parallel_runs=5, hyperparameters_randomtune = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom Ramp Schedules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CIM-AHC, CIM-CAC solvers contain hyperparameters which could vary over the course of one solver run. To set an arbitrary feedback schedule, the arguments custom_feedback_schedule, custom_pump_schedule, and cac_amplitude_ramp can be specified with a specific anonymous function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No External Field Detected\n",
      "Target Ising Energy: -inf.\n",
      "Best Ising Energy Found: -29.0.\n",
      "Corresponding Spin Configuration: [-1. ... -1.].\n",
      "Time Elapsed: 1.629051923751831.\n",
      "Number of Runs Completed: 5.\n"
     ]
    }
   ],
   "source": [
    "timestep_count = 1000\n",
    "my_feedback_schedule = lambda x: np.sin(x/20)#torch.sin(torch.arange(timestep_count))\n",
    "test = Ising(J).solve(cac_gamma=gamma, num_runs=5, num_timesteps_per_run = timestep_count, custom_feedback_schedule=my_feedback_schedule, hyperparameters_randomtune = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting a Target Energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For benchmarking purposes, the ground-state Ising energy of inputted Ising problems can be specified to check whether the solver has reached it. This is modified via the target_energy argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No External Field Detected\n",
      "Target Ising Energy: -29.0.\n",
      "Best Ising Energy Found: -29.0.\n",
      "Corresponding Spin Configuration: [-1. ... -1.].\n",
      "Time Elapsed: 8.247757911682129.\n",
      "Number of Runs Completed: 5.\n",
      "Success Probability: 1.0.\n"
     ]
    }
   ],
   "source": [
    "timestep_count = 5000\n",
    "ground_state_ising_energy = -29.0\n",
    "my_feedback_schedule = lambda x: np.sin(x/33)\n",
    "test = Ising(J).solve(cac_gamma=gamma, num_runs=5, target_energy=ground_state_ising_energy, num_timesteps_per_run = timestep_count, \n",
    "                      custom_feedback_schedule=my_feedback_schedule, hyperparameters_randomtune = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Calculations for an optical Coherent Ising Machine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To translate this simulation into a realistic estimate of how much time this would take on a real FPGA-powered measurement-feedback Coherent Ising Machine, we multiply the number of matrix vector multiplications in our run by the roundtrip time reported in (https://doi.org/10.1126/science.aah5178)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of discrete steps (roundtrips): 5000\n",
      "Roundtrip time: 1.6e-06 second\n",
      "Total CIM runtime: 0.008 second(s)\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of discrete steps (roundtrips): {timestep_count}\")\n",
    "roundtrip_time = 1.6e-6 #Roundtrip time in seconds, as reported in (https://doi.org/10.1126/science.aah5178)\n",
    "print(f\"Roundtrip time: {roundtrip_time} second\")\n",
    "print(f\"Total CIM runtime: {timestep_count*roundtrip_time} second(s)\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.2 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.2"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e8543a852f2c053fa96e40636049c257ee9dfe0580045d9c365f792282b251d7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
